function [P s coeffs]=specw(x,opts);

%function [P s P95 Pc Pprefit]=specw(x,dt,nw,opts);
%
%INPUT
%
%x     equally-spaced time series
%dt    time spacing
%nw    number of windows (for mtm), nw=1 gives a periodogram
%opts: options
%      1. Prewhitening:
%         opts.p, number of autoregressive coefficients (default 1)
%         opts.q, number number of moving-average coefficients (default 0)
%      2. Postwhitening:
%         opts.smoothing; choose from 'gaussian' (default),
%         'lowess', or others available through smoothdata.m
%      3. opts.plot: 1=yes, 0=no (default).
%
%OUTPUT
%
%P     Squared fourier coefficients (nw=0) or power spectral
%      density (nw>0)
%s     Frequency
%P95   95% confidence interval as a function of s
%Pc    Background continuum estimate
%R     Residual, computed as log(P)-log(Pc).


if ~isfield(opts,'p'), opts.p=1; end;
if ~isfield(opts,'q'), opts.q=0; end;
if ~isfield(opts,'plot'), opts.plot=0; end;    

%Expected relationships based on the number of degrees of freedom
k=2*opts.nw-1;
coeffs.dof=2*k;
coeffs.Vexp=psi(1,coeffs.dof/2);  %expected variance of log(P)-log(Pc)
coeffs.bias=psi(0,k)+log(2)-log(k*2); %expected bias in the mean estimate of log(Pc)

%prewhiten using an ARIMA model
if opts.p>0 | opts.q>0, 
    [xpre dummy]=prewhitenPH(x,opts.p,opts.q);
else,
    xpre=x;
end;
coeffs.p=opts.p;
coeffs.q=opts.q;
%coeffs.v=dummy.v;
[dummy s]=pmtmPH(xpre,opts.dt,opts.nw,0); 
P=log(dummy)-(mean(log(dummy))-coeffs.bias); 
coeffs.varP=var(P);

if nargout>1,
    alpha1=0.05*2*k/length(x);
    coeffs.conf=log(gaminv(1-alpha1,k,1/k)); %95 confidence level relative to the log of the background continuum
    
    alpha2=0.01*2*k/length(x);
    coeffs.confb=log(gaminv(1-alpha2,k,1/k)); %99 confidence level for multiple testing regime (Bonferroni correction)
    
end;
